clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%***********PARAMETERS SPECIFICATION Time Varying Var**************
gbit = 1;
maxgbitInit     =  12000;        % total number of draws
burnInit        =  10000;        % burn in period

maxgbit         =   2500;        % total number of draws
burn            =   500;        % burn in period

jump            =     5;           % collect one each jump draws


horizon = 3;
p = 2;


% tight_Q = 0.001;
% tight_W = 0.001;
% tight_S = 0.01;

post_index = 1;
trs=0;
trsMax=10;
MaxEig=1.0;
MaxRoot=MaxEig+1;
c=1;
L=p;





[A B] = xlsread('USdataA.xls',5);
AA=A(1:end,2)./A(1:end,3)*100;
INFL =(log(AA(2:end,1))-log(AA(1:end-1,1)))*100;
%IR =(log(A(2:end,3))-log(A(1:end-1,3)))*400;
M2 =(log(A(2:end,4))-log(A(1:end-1,4)))*100;
GDP =(log(A(2:end,3))-log(A(1:end-1,3)))*100;
DATE =(1870:2008)';
DATA = [INFL M2 GDP];


%******************************************************************
sample_y  = 1870;
s1 = find(DATE(:,1) == sample_y );
DATE = DATE(s1:end,:);
DATA = DATA(s1:end,:);


start_eval_y  = 1900;   
start_sim = find(DATE(:,1) == start_eval_y )-1;

T =size(DATA,1);
%EstWind = 80;
% Z = [INFL GDP];

INFL = DATA(:,1);

EstWind = 10;

tic
count = 0;
for j = start_sim:T-horizon
    j
    count = count +1;
    disp('--------------------------')
    disp('--------------------------')
    disp(['Now running the forecast at ',num2str(DATE(j,1))])




 infl = DATA(1:j,1);
 %gdp = DATA(1:j,3);
m2 = DATA(1:j,2);


    X = [infl m2 ];
%     anto(:,:,j)=cov(X(j-20+1:j,:));
%    aa=cov(X(j-20+1:j,:));
 
 thresh=24;

qq= kron(ones(p,1),kron(std(X)',ones(p,1)));
qq1 = [std(X(:,1)); qq(1:4,1); std(X(:,2)); qq(5:end,1)]; 

qq2=qq1*qq1';

for jj=1:size(qq2,1)
    qq2(jj,jj)=qq1(jj,1);
end
   
   tight_Q = qq2.^(-1/0.32);
 
       tight_Q1 = qq2(1:3,1:3).^(-1/0.38);
    

tight_W = 0.006;
tight_S = 0.006;
    
    PRED_var_rec = VAR_fore(X,p,horizon);
    PRED_var_rol = VAR_fore(X(end-EstWind+1:end,:),p,horizon);


    PRED_ar_rec(1,:) = VAR_fore(X(:,1),p,horizon);
    PRED_ar_rol(1,:) = VAR_fore(X(end-EstWind+1:end,1),p,horizon);




    MixingKSC
    %clear PRED_tv_var PB
    if count == 1
        %% girare con init, chiedere parametri come output, var resid very
        %% important since it is going to be the input for the next start
        [PRED_tv_var, PB,prevR, prevy2star,prevB,prevA] = ...
            TvVarEstKSC(X,tight_Q,tight_W,tight_S,gbit,maxgbitInit,burnInit,jump,post_index,trs,trsMax,MaxEig,MaxRoot,c,L,horizon,ksc);

    else

        [PRED_tv_var, PB, prevR, prevy2star,prevB,prevA] = ...
            TvVarInputEstKSC(X,tight_Q,tight_W,tight_S,gbit,maxgbit,burn,jump,post_index,trs,trsMax,MaxEig,MaxRoot,c,L,horizon,prevR,prevy2star,prevB,prevA,ksc);
   
    end;

    for jg=1:size(PB,3)
        MatA=inv(eye(L*size(X,2))-companion(PB(:,end,jg),L,size(X,2),c));
        LM(:,jg)=MatA(1:size(X,2),1:size(X,2))*PB(1:L*size(X,2)+c:end,end,jg);
    end

    LM_d(j,:) = prctile(squeeze(LM)',16);
    LM_m(j,:) = prctile(squeeze(LM)',50);
    LM_u(j,:) = prctile(squeeze(LM)',84);





    clear PRED_tv_ar


    if count == 1
        [PRED_tv_ar, PBUni,prevRUni, prevy2starUni,prevBUni,prevAUni] = ...
            TvVarEstKSC(X(:,1),tight_Q1,tight_W,tight_S,gbit,maxgbitInit,burnInit,jump,post_index,trs,trsMax,MaxEig,MaxRoot,c,L,horizon,ksc);
    else

           [PRED_tv_ar, PBUni, prevRUni, prevy2starUni,prevBUni,prevAUni] = ...
            TvVarInputEstKSC(X(:,1),tight_Q1,tight_W,tight_S,gbit,maxgbit,burn,jump,post_index,trs,trsMax,MaxEig,MaxRoot,c,L,horizon,prevRUni,prevy2starUni,prevBUni,prevAUni,ksc);
           
    end;

    for jg=1:size(PBUni,3)
        LMUni(:,jg)=(1-ones(1,L)*PBUni(c+1:end,end,jg))^(-1)*PBUni(1,end,jg);
    end

    LM_uni_d(j,1) = prctile(squeeze(LMUni)',16);
    LM_uni_m(j,1) = prctile(squeeze(LMUni)',50);
    LM_uni_u(j,1) = prctile(squeeze(LMUni)',84);

    
    for h = 1:horizon

        fore_var_rec_cum{h}(j+h,:)    = mean(PRED_var_rec(:,1:h),2)';
        fore_var_rol_cum{h}(j+h,:)    = mean(PRED_var_rol(:,1:h),2)';

        fore_ar_rec_cum{h}(j+h,:)    = mean(PRED_ar_rec(:,1:h),2)';
        fore_ar_rol_cum{h}(j+h,:)    = mean(PRED_ar_rol(:,1:h),2)';

        fore_var_tv_d_cum{h}(j+h,:) = prctile(squeeze(mean(PRED_tv_var(:,1:h,:),2))',16);
        fore_var_tv_m_cum{h}(j+h,:) = prctile(squeeze(mean(PRED_tv_var(:,1:h,:),2))',50);
        fore_var_tv_u_cum{h}(j+h,:) = prctile(squeeze(mean(PRED_tv_var(:,1:h,:),2))',84);

        fore_ar_tv_d_cum{h}(j+h,:) = prctile(squeeze(mean(PRED_tv_ar(:,1:h,:),2))',16);
        fore_ar_tv_m_cum{h}(j+h,:) = prctile(squeeze(mean(PRED_tv_ar(:,1:h,:),2))',50);
        fore_ar_tv_u_cum{h}(j+h,:) = prctile(squeeze(mean(PRED_tv_ar(:,1:h,:),2))',84);

        fore_var_tv_m_cum_gibb{h} = squeeze(mean(PRED_tv_var(:,1:h,:),2));
        fore_ar_tv_m_cum_gibb{h} = squeeze(mean(PRED_tv_ar(:,1:h,:),2));
        
        fore_rw{h}(j+h,1)= mean(X(end:end,1));
        true_cum{h}(j+h,1) = mean(INFL(j+1:j+h,1));


    end

      eval(['save ', 'results_var_infl_m2_and_ar_',num2str(j), ' fore_var_tv_m_cum_gibb fore_ar_tv_m_cum_gibb fore_rw true_cum fore_var_rec_cum fore_var_rol_cum fore_ar_rec_cum fore_ar_rol_cum LM_m LM_uni_m LM LMUni' ])



end

